Quantum Monte Carlo calculations of structural properties of FeO under pressure 
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We determine the equation of state of stoichiometric FeO employing the diffusion Monte Carlo 
method. The fermionic nodes are fixed to those of a wave function having the form of a single 
Slater determinant. The calculated ambient pressure properties (lattice constant, bulk modulus and 
cohesive energy) agree very well with available experimental data. At approximately 65 GPa, the 
lattice structure is found to change from rocksalt type (Bf) to NiAs based (inverse B8). 

PACS numbers: 72.80.Ga, 71.15.-m, 71.20.-b, 64.30.+t 



Transition metal oxides are solids with strong electron- 
electron correlations that lead to a rich variety of ob- 
served structural and electronic phases. Almost all tran- 
sition metal oxides are a problem on its own due to 
competitive interplay among correlation and exchange in 
d— subshells, crystal fields, d—p hybridization and charge 
transfer. Elucidation of high-pressure properties of FeO 
is of particular interest in geophysics, since this com- 
pound belongs to constituents of the deep Earth's in- 
terior. FeO is also one of the challenging simple oxides 
due to the nominally open-shell occupation of the 3c? lev- 
els. Indeed, this system proved to be problematic for 
the density functional theory (DFT) in its local density 
or generalized gradient approximations (LDA or GGA). 
For example, both LDA and GGAspredict an incorrect 
ground state lattice structure P, @, 0] . 

At ambient conditions, FeO crystallizes in Bl (NaCl- 
type) structure. It is antiferromagnetically ordered at 
temperatures below 198 K and this ordering is accom- 
panied by a small rhombohedral distortion denoted as 
rBl — the unit cell is stretched along the [111] body diag- 
onal. In shock-wave studies it was observed that around 
70 GPa the oxide transforms to a different structure 
which was inferred as B2 (CsCl-type) in analogy with 
similar materials, but LDA calculations hinted that much 
larger pressure, around 500 GPa, would be needed to 
stabilize B2 against the Bl phase Besides that, no 
such structural transition was detected in static compres- 
sion experiments [|| , unless the material was significantly 
heated up X-ray diffraction performed along the 

high-temperature static compression revealed that the 
high-pressure structure is actually B8 (NiAs-type) 

There are two distinct ways of putting FeO on NiAs 
lattice, the so-called normal B8, where iron occupies Ni 
sites, and inverse B8 (iB8 for short), where iron sits on As 
sites. It is the latter configuration that comes from band- 
structure theories as the more stable of the two Q, Q> HI • 
Also, reinterpretation of the data of Ref. suggests that 
the high pressure phase is a mixture of B8 and iB8 phases 
[l[. On the theoretical side, the introduction of the iB8 
structure into the picture revealed a serious deficiency in 
the LDA (and GGA) as applied to FeO, since the iB8 
phase is predicted more stable not only than B8 but also 



than Bl at all pressures, which contradicts experimental 
findings. It was demonstrated that inclusion of Coulomb 
U to better account for electron-electron correlations al- 
leviates this problem 0, 0] . 

In this Letter, we calculate the equation of state of 
stoichiometric FeO using the fixed- node diffusion Monte 
Carlo method (DMC) [8|, a many-body computational 
technique that accurately treats even strongly correlated 
systems. Based on the aforementioned studies, we con- 
fined ourselves to only two structures — Bl with the type- 
II antiferromagnetic (AFM) ordering (symmetry group 
R3m) and iB8 also in the AFM state (group P6m2). We 
show that DMC, even in its simplest version based on 
a single-determinant Slater-Jastrow wave function, pro- 
vides a very consistent picture of this complicated system 
that closely follows experimental data including estima- 
tion of the transition pressure. 

The guiding wave function that defines (fixes) the 
fermionic nodes in our DMC simulations is of the Slater- 
Jastrow type 4"g = ex p[>^L where 
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The lower-case indices in Eq. (jlbj) run over electrons, 
while the upper-case index denotes ions. The Jastrow 
correlation factor J contains one- and two-body terms, 
g and /, that have the same form as in Ref. and 
represent 17 variational parameters that were optimized 
within variational Monte Carlo (VMC) framework. The 
single determinant of spinorbitals ip a becomes a prod- 
uct of spin-up and -down determinants of spatial orbitals 
<j)p} after fixing the electron spins, = = N/2, 
while the overall state is a spin-unrestricted antiferro- 
magnet. 

The large energy scale of core electrons poses a diffi- 
culty to the DMC in an analogous way as it does to plane- 
wave based electronic structure techniques. Therefore, 
we replace the atomic cores by norm-conserving pseu- 
dopotentials [Toj ] within the so-called localization approx- 
imation 



111 ] . We utilize only small-core pseudopotentials 
to minimize losses in accuracy as much as possible. We 
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have argued recently [12| that even small-core pseudo-po- 
tentials could lead to imprecisions in description of tran- 
sition metal compounds if a spin-related transition oc- 
curs as a part of the phenomena of interest. However, 
it should not affect the present calculations, since iron 
atoms in FeO stay in a high-spin state in the whole range 
of pressures we study in both Bl and iB8 phases. This 
applies to the DMC results as well as to the flavors of 
DFT that we used to construct the Slater determinants. 

The quality of the fixed-node DMC total energy is 
determined solely by the quality of fermionic nodes of 
the guiding wave function. When the form of Eq. |T|) is 
adopted, the parameters controlling location of fermionic 
nodes are the one-electron orbitals {4> a ,<j)p}. The direct 
VMC optimization of these orbitals is currently impracti- 
cal for simulation sizes required in the present study. In- 
stead, we use one-electron orbitals from spin-unrestricted 
calculations with hybrid PBEO functional given as 
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Here E PBE and E PBE are exchange and correlation 
parts of the PBE-GGA E^ F is the exact exchange 
from Hartree-Fock (HF) theory and the weight a is in 
the range < a < 1. We have found that the inclusion 
of exact exchange term into PBE-GGA has similar effect 
as Coulomb U in the LDA+Z7 method. It opens a gap in 
the electronic spectrum of the AFM Bl phase and stabi- 
lizes it relative to the iB8 structure. Both the gap and 
the transition pressure increase with increasing a. With 
mixing weight a — 0.05, the iB8 is still more stable than 
Bl everywhere, the Bl to iB8 transition occurs at 5 GPa 
for a = 0.1 and at 43 GPa for a = 0.2. Note that in MnO, 
which exhibits similar structural transition, the experi- 
mental range of transition pressures is reached already 
for a w 0.1 [IH, while in the present case of FeO, the 
transition takes place much sooner than in experiments 
(> 70 GPa) even for twice as much exact-exchange con- 
tent in PBEO. 

The exchange-correlation functional of Eq. (|2|) 
defines one-parametric class of single-particle or- 
bitals {<j^a\ 4>g } , which can be used to minimize the 
DMC fixed-node error by varying the exact exchange 
weight a. Although in simple insulators, such as silicon, 
the differences between fixed-node energies correspond- 
ing to various sets of one-particle states were found to be 
rather marginal [l5T |. more pronounced differences have 
been obtained for transition metal compounds. In iso- 
lated molecules of transition metal monoxides, the fixed- 
node DMC energies with orbitals from B3LYP (hybrid 
functional similar to PBEO) are noticeably lower than 
with HF or pure DFT orbitals [1, HI]. DMC optimiza- 
tion of the exact exchange proportion in the B3LYP was 
performed in Ref. for the MnO molecule. The opti- 
mal value was reported as approximately 17%, but the 
minimum was rather broad and shallow and values be- 
tween 5% and 30% were almost equivalent. Therefore, 



we have chosen the weight a in PBEO to be a = 0.2, 
corresponding to 20% of GGA exchange being replaced 
with the exact exchange. In the following we abbrevi- 
ate this functional as PBE020- This choice is compatible 
with findings of Ref. 16 and leads to reasonable ambient- 
pressure properties of FeO already within (hybrid) DFT, 
i.e., Bl structure is insulating and more stable than iB8. 
We also checked that at equilibrium the PBEO20 orbitals 
provide DMC energy 0.3 eV per FeO lower than orbitals 
from the HF approximation. 

In our simulations, the infinite crystal was modelled by 
a periodically repeated simulation cell containing 8 FeO 
units, i.e., 176 valence and semi-core electrons. Although 
such a system is certainly not small to deal with in an 
explicitly many-body fashion, it turns out that finite-size 
effects are significant if not treated properly. The origin 
of these finite-size errors is twofold. One part is related to 
incorrect momentum quantization due to confinement of 
electrons into the simulation cell, the second part comes 
from the artificial periodicity of exchange-correlation hole 
due to periodic extension of Coulomb potential using the 
Ewald summation. 

The problem associated with confinement appears also 
in mean-field band theories, where it is exactly resolved 
by integration over the first Brillouin zone, while work- 
ing only within the smallest "simulation" cell possible, 
the primitive cell. Each k-point in the first Brillouin zone 
corresponds to a different boundary condition imposed on 
the primitive cell. The momentum integration is equiv- 
alent to averaging over these so-called twisted boundary 
conditions. Analogous averaging procedure performed in 
a many-body simulation does not represent a complete 
correction, but it pro ved to be very efficient within Monte 
Carlo algorithms [171 ] . In this study we deal only with in- 
sulating states, which simplifies matters considerably and 
we have found that just 8 twists are enough for our simu- 
lation cell size. We have verified that within PBEO20 the 
total energy obtained in our simulation cell using just 8 
k-points differs only ks 0.01 eV per FeO from fully con- 
verged integral over Brillouin zone. 

The second part of finite-size errors, originating from 
artificial periodicity of the Ewald potential, is accounted 



for with the aid of the correction introduced in Ref. 118 . 
The estimate for the energy at infinite volume is written 
as 
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The static structure factor is defined as S(k) = 
(^ol/SkP-kl^o)/^ with pk standing for a Fourier com- 
ponent of the electron density. The integral in Eq. ([3]) 
runs over a domain D centered around k = and hav- 
ing volume 87r 3 /f2, where ft is volume of the simulation 
cell. The structure factor S(k) is evaluated within the 
DMC at a discrete set of points and then extrapolated 
towards k = 0. Performance of the correction given by 
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FIG. 1: (color online) Finite-size errors of the twist averaged 
DMC energy at volumes 15.9 A 3 /FeO (left) and 20.4 A 3 /FeO 
(right). Pure Ewald energies are shown as red squares, values 
corrected according to Eq. ([3| are represented by blue circles. 
Note that the finite-size errors increase with increasing the 
electron density. The energy in infinite cell Eoo is extrapo- 
lated from the data shown. Statistical errorbars are smaller 
than symbol sizes except for the largest cell in the left panel. 



Eq. © applied to FeO is illustrated in Fig. [TJ where we 
plot the total energy at two different electron densities 
as calculated in simulation cells of varied size up to 16 
FeO units, i.e., 352 valence and semi-core electrons. The 
correction removes more than 90% of the finite-size error 
introduced by the periodic electron-electron interaction 
potential and enables us to replace expensive size scaling 
analysis of Fig. Q] by a simple formula, Eq. ([3]). 

Results. First we discuss properties of the Bl phase 
around equilibrium volume. For the cohesion energy, 
E C oh = E-p c + Eq — Ef c o, our DMC simulations yield 
9.66 ± 0.04 eV/FeO that matches 9.7 eV/FeO deduced 
from experimental formation enthalpies [1 91 ] . The elec- 
tronic gap, which we calculate as a difference between 
total energies of the ground state and the first excited 
state at the T— point in our 8 FeO simulation cell, comes 
out as 2.8 ± 0.3 eV. This value is not too far from opti- 
cal absorption edge observed near 2.4 eV [20]. The weak 
feature displayed between 1.0 and 1.5 eV in these experi- 
ments is not reproduced in the picture of FeO we present 
here. However, it is quite possibly related to imperfec- 
tions in structural or magnetic order, since essentially the 
same absorption band was repeatedly observed in other 
systems where Fe atoms act as impurities 
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The DMC energy as a function of volume, together 
with fitted Murnaghan equation of state, is presented in 
Fig. [2j The parameters of the least-square fit are com- 
pared with predictions of other electronic structure meth- 
ods and with experiments in Tab. HI The DMC estimate 
for equilibrium lattice constant ao is in excellent agree- 
ment with experimental value extrapolated to stoichio- 
metric FeO and offers more consistent prediction than 
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FIG. 2: (color online) Total energies of Bl (red squares) and 
iB8 (blue circles) phases. Statistical errorbars are smaller 
than symbol sizes. Lines are least-square fits with Murnaghan 
equation of state. The dashed double tangent corresponds to 
the transition pressure of 65 GPa. Inset shows the data for 
Bl phase over wider volume region, including equilibrium. 



LDA or GGA. On the other hand, the hybrid PBE0 20 
functional, which we used to construct the DMC guid- 
ing wave function, provides a similar value. All methods 
shown in Tab. U provide essentially the same value of bulk 
modulus, which is noticeably larger than typical exper- 
imental reports. The extrapolation to stoichiometry is, 
however, expected to lead to values in the vicinity of the 
theoretical data 

[IHIII. The isothermal pressure deriva- 
tive of the bulk modulus, K' , turns out to be larger in 



DMC than in the density-functional approaches, which 
makes it compatible with elastic- wave experiments [24j . 

The equation of state calculated within diffusion Monte 
Carlo up to large hydrostatic pressures is shown in Fig. [2] 



TABLE I: Equilibrium lattice constant ao, bulk modulus Kq 
and its derivative K' = (dKo/dP)r calculated in this work 
(DMC and PBEO20) compared to selected theories and exper- 
iments. The experimental lattice constant is extrapolated to 
the stoichiometric FeO, whereas the values of Ko and K' are 
not. The extrapolated value of bulk modulus Ko is estimated 
around 180 GPa 





a () (A) 


Ko (GPa) 


K 


DMC 


4.324(6) 


170(10) 


5.3(7) 


PBEO20 


4.328 


182 


3.7 


GGA [3] 


4.28 


180 


3.6 


LDA [5] 


4.136 


173 


4.2 


experiment 


4.334 [23] 


152.3 [24] 


4.92 [24] 
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FIG. 3: (color online) P(V) curves of Bl (red) and iB8 (blue) 
phases. Lines are our DMC fits as in Fig. [2] Points are ex- 
perimental data: filled circles [(| (BJL Feo.gsO), empty circles 
[H (Bl, Fe .94O), empty squares 0] (iB8, Feo.gsO). All Bl 
data [1, 28] were taken at room temperature, the iB8 data 
[3 correspond to 900 K. Bl datasets are shown relative to 
the equilibrium volumes reported in the individual studies to 
approximately remove the non-stoichiometry effects. 



The c/a ratio in the hexagonal iB8 phase, stable at high 
pressures, was optimized within PBEO20 It was found 
to increase from 1.93 at volume 17 A 3 /FeO to 2.03 at 
14 A 3 /FeO. These ratios agree well with experimental 
c/a = 2.01 at 14.83 A 3 /FeO reported in Ref.0. 

The Bl phase, stable at low pressures, was assumed cu- 
bic, i.e., we neglected the rhombohedral distortion. We 
did not use the DFT optimized geometries in this case, 
because DFT based techniques are not conclusive in de- 
tcrmination of the rhombohedral distortion 
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We checked the impact of fixing the cubic symmetry by 
comparing DMC total energies for different distortions 
at high compression, where the distortion has the largest 
impact. At the volume of 15.3 A 3 /FeO, the energy gain 
associated with the rhombohedral distortion was of the 
order of statistical errorbars « 0.02 eV/FeO, i.e., too 
small to affect the results. 

The critical pressure P c of the structural transition 
from Bl to iB8 phase was determined from equality of 
Gibbs potentials, Gbi(P c ) — GiBs(Pc), graphical equiv- 
alent of which — the double tangent — is shown in Fig. [2] 
The value is P c = 65 ± 5 GPa with the errorbar given by 
statistical fluctuations of the Monte Carlo simulations. 
The prediction P c = 65 GPa agrees quite well with shock- 
wave data and with high-temperature static compression 
experiments, except for the fact that our investigation is 
performed at T = 0, for which experiments suggest con- 
siderably higher transition pressure. Our finding could 
be interpreted as an indirect support for existence of an 
energy barrier between the two phases that requires a 
thermal activation to be overcome. 



Another means of comparison with experiments is 
looking at the P(V) equation of state, Fig. [3l Agree- 
ment between data for Bl structure (after extrapolation 
to stoichiometric FeO) is very good, which is in concord 
with similarly good correspondence of ambient pressure 
parameters compared in Tab. [I] Our curve for iB8 struc- 
ture also follows the x-ray data of Ref. rather nicely 
(no stoichiometry related correction was attempted in 
this case). Experimentally, this phase is perhaps some- 
what stiffer than in our calculations, which signals that 
DMC is likely to slightly underestimate the transition 
pressure P c . 

In summary, the equation of state and basic electronic 
structure of FeO calculated with the diffusion Monte 
Carlo agrees very well with many aspects of available 
experimental data. Considering that the DMC is essen- 
tially a parameter-free method and that the fixed-node 
condition was enforced with the aid of a very simple wave 
function (single Slater determinant), the degree of con- 
sistency of the provided picture is quite remarkable. 

We acknowledge support by NSF DMR-0121361 and 
EAR-0530110 grants. This study was enabled by IN- 
CITE allocation at ORNL and by allocation at NCSA. 
QMC simulations were done using QWalk code [29j |. 
and the o ne-p article orbitals were calculated with Crys- 
TAL2003 0. 
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